Ao final desta aula o aluno deve ser capaz de:
Suponha que foi feito um levantamento sobre o número de drogas em uso por gestantes .
DrgGrv <- read.table("Drogas_Gravidez.txt", header=TRUE, sep="\t")
print(DrgGrv)## Drogas Pacientes
## 1 0 1425
## 2 1 1351
## 3 2 793
## 4 3 348
## 5 4 156
## 6 5 58
## 7 6 28
## 8 7 15
## 9 8 6
## 10 9 3
## 11 10 1
## 12 11 0
## 13 12 1
O arquivo Drogas_Gravidez.txt tem duas colunas:
A variável DrgGrv é um data frame:
is.data.frame(DrgGrv)## [1] TRUE
plot (DrgGrv)plot (DrgGrv,type="h")Como foi que eu descobri isto? Quick-R: Line charts
plot (DrgGrv,type="h")
points(DrgGrv)Encontrei mais em Quick-R: Graphical Parameters
plot (DrgGrv,type="h", col="#994F88", lwd=3)
points(DrgGrv, pch=21, col="#994F88", bg="#994F88")source("friendlycolor.R")
plot (DrgGrv,type="h", col=friendlycolor(8), lwd=3)
points(DrgGrv, pch=21, col=friendlycolor(8), bg=friendlycolor(12))Exibindo todas as cores da função friendlycolor():
# friendlydemo.R
# exibe as 46 cores disponíveis em friendlycolor.R
source("friendlycolor.R")
indice <- 1 # primeiro indice
# plota um gráfico vazio
plot(NA, xlim=c(0,9), ylim=c(0,7),
xlab=NA, ylab=NA, axes = FALSE)
for (x in 1:8)
{
for (y in 1:6)
{
# pula se acabaram as cores
if (indice > 46) {next}
points ( x, 7-y,
pch=21, cex=5,
col="black",
bg=friendlycolor(indice)
)
cortexto <- "black"
if (
(indice >= 25 & indice <= 27)
|
(indice >= 31 & indice <= 36)
)
{
cortexto <- "white"
}
text(x, 7-y, indice, col=cortexto)
# incrementa (proxima cor)
indice <- indice+1
}
}DrgGrv <- read.table("Drogas_Gravidez.txt", header=TRUE, sep="\t")
print(DrgGrv)## Drogas Pacientes
## 1 0 1425
## 2 1 1351
## 3 2 793
## 4 3 348
## 5 4 156
## 6 5 58
## 7 6 28
## 8 7 15
## 9 8 6
## 10 9 3
## 11 10 1
## 12 11 0
## 13 12 1
names(DrgGrv)## [1] "Drogas" "Pacientes"
DrgGrv$Porcentagem <- round(DrgGrv$Pacientes/sum(DrgGrv$Pacientes)*100,2)
names(DrgGrv)## [1] "Drogas" "Pacientes" "Porcentagem"
print(DrgGrv)## Drogas Pacientes Porcentagem
## 1 0 1425 34.05
## 2 1 1351 32.28
## 3 2 793 18.95
## 4 3 348 8.32
## 5 4 156 3.73
## 6 5 58 1.39
## 7 6 28 0.67
## 8 7 15 0.36
## 9 8 6 0.14
## 10 9 3 0.07
## 11 10 1 0.02
## 12 11 0 0.00
## 13 12 1 0.02
sum(DrgGrv$Porcentagem)## [1] 100
DrgGrv$Porcentagem < 0## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE
sum(DrgGrv$Porcentagem < 0)## [1] 0
DrgGrv$Porcentagem >= 0 & DrgGrv$Porcentagem <= 100## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
sum (DrgGrv$Porcentagem >= 0 & DrgGrv$Porcentagem <= 100)## [1] 13
Como DrgGrv agora tem três colunas, precisamos indicar quais colunas quero usar como \(x\) (DrgGrv$Drogas) e \(y\) (DrgGrv$Porcentagem)
source("friendlycolor.R")
plot (DrgGrv$Drogas, DrgGrv$Porcentagem,
type="h", col=friendlycolor(20), lwd=3)
points(DrgGrv$Drogas, DrgGrv$Porcentagem,
pch=21,
col=friendlycolor(20), bg=friendlycolor(23))… consequentemente, os nomes das variáveis aparecem nos eixos. Isto é corrigido assim:
source("friendlycolor.R")
plot (DrgGrv$Drogas, DrgGrv$Porcentagem,
type="h", col=friendlycolor(20), lwd=3,
xlab="Número de drogas utilizadas",
ylab="Porcentagem de gestantes")
points(DrgGrv$Drogas, DrgGrv$Porcentagem,
pch=21,
col=friendlycolor(20), bg=friendlycolor(23))O número de drogas utilizadas é uma variável quantitativa discreta.
Apenas para demonstrar (não é a forma apropriada), caso a variável fosse quantitativa contínua, poderíamos usar um density plot
Primeiro criamos uma variável para conter a densidade de probabilidades
# desdobra a quantidade de pacientes
populacao <- c() # cria um vetor vazio
for (r in 1:nrow(DrgGrv))
{
# acumula no vetor o número de drogas
populacao <- c(populacao, rep(DrgGrv$Drogas, times=DrgGrv$Pacientes))
}
# usa a funcao R que transforma em densidade de probabilidade
densidade <- density(populacao)source("friendlycolor.R")
plot (densidade,
col=friendlycolor(20), lwd=3,
xlab="Número de drogas utilizadas",
ylab="Porcentagem de gestantes")A função plot() tem muitos parâmetros; aqui usamos o parâmetro main.
source("friendlycolor.R")
# note o titulo em duas linhas usando-se \\n
plot (densidade,
main = "Distribuição do uso de\ndrogas em gestantes",
col=friendlycolor(20), lwd=3,
xlab="Número de drogas utilizadas",
ylab="Porcentagem de gestantes")Tudo o que foi feito até agora pode ser colocado em um único Rscript
# gestantes.R
# le os dados e cria os graficos
source("friendlycolor.R")
# le os dados
DrgGrv <- read.table("Drogas_Gravidez.txt",
header=TRUE, sep="\t")
# cria coluna para a porcentagem
DrgGrv$Porcentagem <- round(DrgGrv$Pacientes/
sum(DrgGrv$Pacientes)*100,2)
# exibe a tabela de dados
cat("Utilização de drogas em gestantes\n")
print(DrgGrv)
# exibe o grafico no estilo 'histograma'
# (variavel quantitativa discreta)
plot (DrgGrv$Drogas, DrgGrv$Porcentagem,
main = "Distribuição do uso de\ndrogas em gestantes",
type="h", col=friendlycolor(20), lwd=3,
xlab="Número de drogas utilizadas",
ylab="Porcentagem de gestantes")
points(DrgGrv$Drogas, DrgGrv$Porcentagem,
pch=21,
col=friendlycolor(20), bg=friendlycolor(23))
# tratando a variavel como quantitativa continua
# (apenas por exemplo; nao eh apropriado fazer isto)
# cria uma populacao com as quantidades
# de drogas utilizadas pelas pacientes
populacao <- c() # cria um vetor vazio
for (r in 1:nrow(DrgGrv))
{
# acumula no vetor o número de drogas
populacao <- c(populacao,
rep(DrgGrv$Drogas, times=DrgGrv$Pacientes))
}
# usa a funcao R que transforma
# em densidade de probabilidade
densidade <- density(populacao)
# exibe o grafico no estilo 'density plot'
plot (densidade,
main = "Distribuição do uso de\ndrogas em gestantes",
col=friendlycolor(20), lwd=3,
xlab="Número de drogas utilizadas",
ylab="Porcentagem de gestantes")## Utilização de drogas em gestantes
## Drogas Pacientes Porcentagem
## 1 0 1425 34.05
## 2 1 1351 32.28
## 3 2 793 18.95
## 4 3 348 8.32
## 5 4 156 3.73
## 6 5 58 1.39
## 7 6 28 0.67
## 8 7 15 0.36
## 9 8 6 0.14
## 10 9 3 0.07
## 11 10 1 0.02
## 12 11 0 0.00
## 13 12 1 0.02
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
|
A função R é dbinom(x, size, prob), indicando, respectivamente, quantas jogadas, o total de jogadas e a probabilidade de sucesso de uma jogada. Para uma moeda balanceada (prob=0.5), a probabilidade de 0 sucesso (x=0) em 3 jogadas (size=3) é:
1 sucesso em 3 jogadas:
2 sucesso em 3 jogadas:
3 sucesso em 3 jogadas:
|
Os gráficos foram produzidos utilizando o seguinte código:
source("friendlycolor.R")
jogadas <- 5
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,0.5)
plot(sucesso, probabilidade,
main = paste("Binomial: ",
jogadas, " jogadas", sep=""),
ylim = c(0,0.5),
type="h",
col=friendlycolor(8), lwd=3)
points(sucesso,probabilidade, pch=21,
col=friendlycolor(8),
bg=friendlycolor(12))|
É possível ver todos os valores em uma tabela:
… ou, mais facilmente ainda, criando um data frame:
… as colunas do data frame podem ser renomeadas:
|
Este código foi modificado para usar um data frame.
source("friendlycolor.R")
jogadas <- 15
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,0.5)
binomial <- data.frame(sucesso,probabilidade)
names(binomial) <- c("Sucesso","FR")
plot(binomial$Sucesso,
binomial$FR,
main = paste("Binomial: ",
jogadas, " jogadas", sep=""),
ylim = c(0,0.5),
type="h",
col=friendlycolor(8), lwd=3)
points(binomial$Sucesso,
binomial$FR,
pch=21,
col=friendlycolor(8),
bg=friendlycolor(12))| Observe que a soma de todas as colunas é igual a 1. Então, quanto mais jogadas, maior a dispesão e menor a altura das distribuições. |
|
|
plot(binomial$Sucesso,
binomial$FR,
main = paste("Binomial: ",jogadas, " jogadas", sep=""),
ylim = c(0,max(binomial$FR)),
type="h",
col=friendlycolor(8), lwd=3)
points(binomial$Sucesso,
binomial$FR,
pch=21,
col=friendlycolor(8),
bg=friendlycolor(12))source("friendlycolor.R")
p.sucesso <- 0.7 # *** probabilidade de sucesso ***
jogadas <- 15
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,p.sucesso)
plot(sucesso, probabilidade,
main = paste("Binomial: ",
jogadas, " jogadas, ",
"P[s] = ", p.sucesso,
sep=""),
ylim = c(0,max(probabilidade)),
type="h",
col=friendlycolor(20), lwd=3)
points(sucesso,probabilidade, pch=21,
col=friendlycolor(20),
bg=friendlycolor(23))source("friendlycolor.R")
p.sucesso <- 0.5 # *** probabilidade de sucesso ***
jogadas <- 15
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,p.sucesso)
# criando um data frame
binomial <- data.frame(sucesso,probabilidade)
names(binomial) <- c("Sucesso", "FR")
# grafico
cauda <- 1 # a partir de onde hachurar
# largura da hachura
hachura <- 500/jogadas
if (hachura < 10) {hachura <- 10}
if (hachura > 35) {hachura <- 35}
plot(binomial$Sucesso,
binomial$FR, # *** usando o dataframe
main = paste("Binomial: ",
jogadas, " jogadas, ",
"P[sucesso] = ", p.sucesso,
sep=""),
ylim = c(0,max(probabilidade)),
type="h",
col=friendlycolor(8), lwd=3)
points(sucesso,probabilidade, pch=21,
col=friendlycolor(8),
bg=friendlycolor(12))
# cauda esquerda
lines (binomial$Sucesso[binomial$Sucesso<cauda],
binomial$FR[binomial$Sucesso<cauda],
col=paste(friendlycolor(8),"88",sep=""),
lwd=hachura, type="h")
# cauda direita
lines (binomial$Sucesso[binomial$Sucesso>15-cauda],
binomial$FR[binomial$Sucesso>15-cauda],
col=paste(friendlycolor(8),"88",sep=""),
lwd=hachura, type="h")total <- sum(dbinom(c(0,15),jogadas,p.sucesso))
cat("Total = ", total, sep="")## Total = 6.103516e-05
\(P[s \le 0] + P[s \ge 15] \approx 6/10000\)
|
Note que o data frame facilitou a construção do gráfico e hachura de suas caudas. No entanto, os eixos dos gráficos receberam o nome das variáveis (o que não é conveniente). Dois parâmetros de plot() resolvem o problema: |
total <- sum(dbinom(c(0:1,14:15),jogadas,p.sucesso))
cat("Total = ", total, sep="")## Total = 0.0009765625
\(P[s \le 1] + P[s \ge 14] \approx 1/1000\)
total <- sum(dbinom(c(0:2,13:15),jogadas,p.sucesso))
cat("Total = ", total, sep="")## Total = 0.007385254
\(P[s \le 2] + P[s \ge 13] \approx 0.7\%\)
total <- sum(dbinom(c(0:3,12:15),jogadas,p.sucesso))
cat("Total = ", total, sep="")## Total = 0.03515625
\(P[s \le 3] + P[s \ge 12] \approx 3.52\%\)
total <- sum(dbinom(c(0:4,11:15),jogadas,p.sucesso))
cat("Total = ", total, sep="")## Total = 0.1184692
\(P[s \le 4] + P[s \ge 11] \approx 11.85\%\)
total <- sum(dbinom(c(0:5,10:15),jogadas,p.sucesso))
cat("Total = ", total, sep="")## Total = 0.3017578
\(P[s \le 5] + P[s \ge 10] \approx 30.18\%\)
\(P[s \le 3] + P[s \ge 12] \approx 3.52\%\)
\(H_0: P[sucesso] = 0.5\)
\(H_1: P[sucesso] \ne 0.5\)
|
\(\alpha\) … probabilidade do erro do tipo I (rejeitar \(H_0\) incorretamente) |
Dado um valor a receber em moedas de R$1.00, metade da quantia eh oferecida em moedas com um balanceamento de referencia, e metade em moedas falsas, conhecidas por balanceamento distinto. Seu desafio eh distinguir os dois conjuntos atraves de experimentos. Numero de moedas (inteiro, default=10000): 10000
Para testar se a moeda eh verdadeira, joga-se cara ou coroa certo numero de vezes cada moeda (um experimento). Numero de lancamentos por experimento (numero inteiro, default=15): 15
Qual a proporcao maxima de moedas verdadeiras que voce aceita perder, i.e. alfa = probabilidade do erro do tipo I ou de falso-positivo). (numero entre 0 e 1). alfa (default=0.05): 0.05
As moedas verdadeiras tem balanceamento de referencia (H0). (caso queira moedas balanceadas, escolha o valor igual a 0.5) Qual a probabilidade de sortear coroa para uma moeda verdadeira? (número entre 0 e 1). P[coroa|H0] (default=0.5): 0.5
source("Goodcoin.R")Certo método educacional consegue ensinar higiene pessoal a 50% dos pacientes com autismo.
Com o novo método proposto pelo Instituto Ayres Soares, no entanto, entre 15 crianças acompanhadas, 10 (66.66%) conseguiram aprender a cuidar de sua higiene.
Se os métodos fossem iguais, somente metade das crianças (7 ou 8) deveriam aprender.
O novo método é melhor?
\(H_0: \mu_{novo} = 0.5\)
\(H_1: \mu_{novo} > 0.5\)
\(\alpha = 0.05 = 5\%\)
(este teste é unicaudal, tem direção)
Probabilidade de observar a melhora de 10 crianças em 15 testadas, sob \(H_0\) (i.e., assumindo-se que o novo tratamento tem o mesmo efeito que o tratamento antigo).
p <- sum(binomial$FR[binomial$Sucesso>=10])
cat("p = ",p,"\n")## p = 0.1508789
\(P[s \ge 10] \approx 15.09\%\)
Probabilidade do erro do tipo I: critério escolhido pelo pesquisador.
alfa <- sum(binomial$FR[binomial$Sucesso>=12])
cat("alfa = ",alfa,"\n")## alfa = 0.01757812
\(P[s \le 3] + P[s \ge 12] \approx 1.78\%\)
Decisão: não se rejeita \(H_0\)
Portanto, não há evidência de que o novo tratamento seja superior ao tratamento tradicional, tomado como referência.
source("Goodcoin.R")Dado um valor a receber em moedas de R$1.00, metade da quantia eh oferecida em moedas com um balanceamento de referencia, e metade em moedas falsas, conhecidas por balanceamento distinto. Seu desafio eh distinguir os dois conjuntos atraves de experimentos. Numero de moedas (inteiro, default=10000): 20000
Para testar se a moeda eh verdadeira, joga-se cara ou coroa certo numero de vezes cada moeda (um experimento). Numero de lancamentos por experimento (numero inteiro, default=15): 15
Qual a proporcao maxima de moedas verdadeiras que voce aceita perder, i.e. alfa = probabilidade do erro do tipo I ou de falso-positivo). (numero entre 0 e 1). alfa (default=0.05): 0.05
As moedas verdadeiras tem balanceamento de referencia (H0). (caso queira moedas balanceadas, escolha o valor igual a 0.5) Qual a probabilidade de sortear coroa para uma moeda verdadeira? (número entre 0 e 1). P[coroa|H0] (default=0.5): 0.5
As moedas falsas tem outro balanceamento. (para simular, forneça uma probabilidade diferente de 0.5 ou deixe em branco para simular somente a moeda verdadeira) Qual a probabilidade de sortear coroa para uma moeda falsa? (número entre 0 e 1). P[coroa|H1] (default=0.5): 0.6666
|
\(\beta\) … probabilidade do erro do tipo II (não rejeitar \(H_0\) incorretamente) |
| \(H_0\) verdadeira | \(H_0\) falsa | |
| não rejeita \(H_0\) | ok | \(\beta\) |
| rejeita \(H_0\) | \(\alpha\) | ok |
Probabilidade de rejeitar \(H_0\) corretamente:
|
O efeito não existe |
O efeito existe |
|
|
sem evidência de efeito |
não rejeitou \(H_0\), corretamente |
\(\beta\) |
|
com evidência de efeito |
\(\alpha\) |
rejeitou \(H_0\), corretamente \(poder = 1-\beta\) |
| Não sabemos se o efeito existe | |
|
não rejeitou \(H_0\), sem evidência de efeito, então… |
… o efeito não existe OU … o efeito existe e a probabilidade de decisão errada é \(\beta\) (se o efeito existe -> erro do tipo II) |
|
rejeitou \(H_0\), com evidência de efeito, então… |
… o efeito não existe e a probabilidade de decisão errada é \(\alpha\) (se o efeito não existe -> erro do tipo I) OU … o efeito existe e a probabilidade de decisão correta é \(1-\beta\), i.e., o poder do teste (de declarar a diferença acertadamente). |
Então, neste exemplo, como não rejeitamos \(H_0\), declarar que os dois tratamentos são iguais (i.e., aceitar \(H_0\)) é a decisão incorreta com probabilidade de …
p.sucesso <- 0.6666
jogadas <- 15
beta <- sum(dbinom(0:11,jogadas,p.sucesso))
cat("beta = ",beta,"\n")## beta = 0.7909156
Não temos evidência para dizer que os dois tratamentos são diferentes, e menos ainda podemos afirmar que os dois são iguais.
A literatura costuma usar o nível de significância \(\alpha=0.05=5\%\) e poder entre \(80\%\) e \(90\%\) (de \(\beta=0.2=20\%\) a \(\beta=0.1=10\%\)).
source("Goodcoin.R")Mesmo com \(alfa=0.6\) …
## cutoff = 7
cat("beta = ",round(beta*100,2),"%\n")## beta = 8.82 %
cat("poder = ",round((1-beta)*100,2),"%\n")## poder = 91.18 %
source("Goodcoin.R")Com 78 crianças e \(alfa=0.05\) …
## cutoff = 46
cat("beta = ",round(beta*100,2),"%\n")## beta = 9.46 %
cat("poder = ",round((1-beta)*100,2),"%\n")## poder = 90.54 %
Aplicando-se o novo método aplicado a 78 crianças:
se 47 ou mais crianças aprenderem higiene pessoal, rejeitamos \(H_0\) e podemos afirmar que o novo tratamento é melhor que o tradicional com \(90\%\) de probabilidade de estarmos corretos (no entanto, pode ser que os tratamentos sejam iguais e estejamos errados, com \(\alpha=5\%\) de probabilidade).
se somente até 46 crianças aprenderem higiene pessoal, não podemos rejeitar \(H_0\). No entanto, como escolhemos poder de \(90\%\) antes do iniciarmos o estudo, podemos aceitar \(H_0\) e afirmar que os dois tratamentos são iguais porque a igualdade não deve ser decorrente de insuficiência amostral (no entanto, podemos estarm enganados e, na verdade, o novo tratamento ser melhor, com \(\beta=10\%\) de probabilidade).
|
Notou que \({47 \over 78} = 0.6025641 \approx 60\%\) e antes, quando foi inconclusivo, \({10 \over 15} = 0.6666667 \approx 67\%\) ? |
source("Goodcoin.R")Caso eu não tenha como avaliar mais do que 15 crianças…
jogadas <- 15
sucesso <- 0:jogadas
p.sucesso <- 0.88 # *** um novo tratamento, com efeito maior
probabilidade <- dbinom(sucesso,jogadas,p.sucesso)
binomial <- data.frame(sucesso,probabilidade)
names(binomial) <- c("Sucesso","FR")
binomial$FA <- NA
binomial$FAdec <- NA
for(b in 1:(jogadas+1))
{
binomial$FA[b] <- sum(binomial$FR[binomial$Sucesso<=b-1])
binomial$FAdec[b] <- sum(binomial$FR[binomial$Sucesso>=b-1])
}
print(binomial)## Sucesso FR FA FAdec
## 1 0 1.540702e-14 1.540702e-14 1.0000000
## 2 1 1.694772e-12 1.710179e-12 1.0000000
## 3 2 8.699832e-11 8.870849e-11 1.0000000
## 4 3 2.764613e-09 2.853322e-09 1.0000000
## 5 4 6.082149e-08 6.367481e-08 1.0000000
## 6 5 9.812534e-07 1.044928e-06 0.9999999
## 7 6 1.199310e-05 1.303802e-05 0.9999990
## 8 7 1.130778e-04 1.261158e-04 0.9999870
## 9 8 8.292370e-04 9.553528e-04 0.9998739
## 10 9 4.729722e-03 5.685075e-03 0.9990446
## 11 10 2.081078e-02 2.649585e-02 0.9943149
## 12 11 6.936925e-02 9.586511e-02 0.9735041
## 13 12 1.695693e-01 2.654344e-01 0.9041349
## 14 13 2.869634e-01 5.523978e-01 0.7345656
## 15 14 3.006283e-01 8.530261e-01 0.4476022
## 16 15 1.469739e-01 1.000000e+00 0.1469739
## cutoff = 11
cat("beta = ",round(beta*100,2),"%\n")## beta = 9.59 %
cat("poder = ",round((1-beta)*100,2),"%\n")## poder = 90.41 %
Com apenas 15 crianças disponíveis, o estudo de um novo método poderá ser conclusivo se o método for capaz de ensinar com sucesso \(88\%\) ou mais das crianças (i.e., tenha efeito maior que os 2/3 do exemplo anterior); neste caso poderá ser considerado igual ao tradicional no caso de até 11 crianças ou melhor que o tradicional se 12 ou mais crianças aprenderem higiene pessoal.
Esta é uma distribuição de probabilidades para uma variável quantitativa discreta (contagem de eventos) quando:
Certo estudo demonstrou que a distribuição mensal de acidentes de trabalho entre moradores de determinada comunidade, entre 1977 e 1987, obedecia uma distribuição de Poisson com média 33 ocorrências por ano.
Qual a probabilidade estimada de que num certo mês do ano sejam observados 3 acidentes de trabalho?
|
A distribuição de Poisson tem um único parâmetro, \(\lambda\), que é a taxa de ocorrência dos eventos. Neste exemplo, a taxa mensal é \({33 \over 12} = 2.75\) |
|
A função R é dpois(x, lambda) indicando, respectivamente, quantas ocorrências e a taxa das ocorrências. No caso de ocorrerem 2.75 (lambda=2.75) eventos para o tempo e o tamanho da população considerada, a probabilidade de não ocorrer evento algum (x=0) é:
de 1 ocorrência:
de 2 ocorrências:
etc. Diferentemente da distribuição binomial, a de Poisson pode ser computada ao infinito. |
Neste exemplo, a ocorrência de exatamente 3 eventos é:
dpois(x=3, lambda=2.75)## [1] 0.2215833
A probabilidade de ocorrer até 3 eventos é a soma das probabilidades de ocorrência de nenhum, 1, 2 ou 3 eventos:
sum(dpois(0:3, lambda=2.75))## [1] 0.70304
Como não há limite superior, caso quiséssemos saber a probabilidade de 4 ou mais eventos, precisamos usar o complemento:
1-sum(dpois(0:3, lambda=2.75))## [1] 0.29696
O seguinte código mostra a aparência desta distribuição:
# Poisson_com_caudas.R
# altere o que deseja
lambda <- 2.75
max_plotar <- 15 # valor maximo para plotar (Poisson vai a + infinito)
cauda <- 3 # a partir de onde hachurar e calcular p da area
cor = "#1965B0" # RGB
# grafico e calculo
hachura <- 500/max_plotar
if (hachura < 10) {hachura <- 10}
if (hachura > 35) {hachura <- 35}
cor_transparencia <- paste(cor,"88",sep="")
eventos <- 0:max_plotar
probs <- dpois(eventos,lambda)
# data frame com a distribuicao de Poisson
poisson <- data.frame(eventos,probs)
names(poisson) <- c("Eventos", "Probabilidade")
plot (poisson$Eventos, poisson$Probabilidade,
main="Distribuição de Poisson",
xlab="Número de ocorrências",
ylab="Probabilidade",
lwd=3, col=cor, type = "h")
points(poisson$Eventos, poisson$Probabilidade,
col=cor, bg=cor, pch=21)
# cauda direita
lines (poisson$Eventos[poisson$Eventos>cauda],
poisson$Probabilidade[poisson$Eventos>cauda],
col=paste(cor,"88",sep=""),
lwd=hachura, type="h")print(poisson)## Eventos Probabilidade
## 1 0 6.392786e-02
## 2 1 1.758016e-01
## 3 2 2.417272e-01
## 4 3 2.215833e-01
## 5 4 1.523385e-01
## 6 5 8.378618e-02
## 7 6 3.840200e-02
## 8 7 1.508650e-02
## 9 8 5.185984e-03
## 10 9 1.584606e-03
## 11 10 4.357667e-04
## 12 11 1.089417e-04
## 13 12 2.496580e-05
## 14 13 5.281228e-06
## 15 14 1.037384e-06
## 16 15 1.901871e-07
prob_dir <- sum(poisson$Probabilidade[poisson$Eventos<=cauda])
cat("\nP[eventos <= ",cauda,"] =",prob_dir)##
## P[eventos <= 3 ] = 0.70304
cat("\nP[eventos > ",cauda,"] =",1-prob_dir)##
## P[eventos > 3 ] = 0.29696
cat("\n")| A distribuição normal tem dois parâmetros, média (\(\mu\)) e desvio-padrão (\(\sigma\)). |
Vamos assumir uma distribuição \(N(\mu=15,\sigma=8)\):
É uma distribuição simétrica, portanto média, moda e mediana coincidem:
Metade da área sob a curva está à esquerda e metade à direita da média.
cerca de 68% da área entre -1 e +1 desvio-padrão:
cerca de 95% da área entre -2 e +2 desvio-padrão:
cerca de 99.7% da área entre -3 e +3 desvio-padrão:
Para cada par de parâmetros \(\mu\) e \(\sigma\), define-se completamente uma distribuição normal.
Admita (não tome estes valores como variáveis da prática médica) que as seguintes variáveis tenham distribuições normais:
# GraficoNormal.R
source ("friendlycolor.R")
variavel <- "Creatinina"
unidade <- "mg/dl"
media <- 0.95
desvpad <- 0.125
x <- seq(from=media-5*desvpad, to=media+5*desvpad, by=0.01)
y <- dnorm(x, mean=media, sd=desvpad)
xy <- data.frame(x,y)
plot(x,y,
main=paste("N(",media,",",desvpad,")",sep=""),
xlab=paste(variavel," (",unidade,")",sep=""),
ylab=NA,
xlim=c(media-4*desvpad,media+4*desvpad),
col=friendlycolor(2),type="l",lwd=2
)… com o mesmo código R, geramos:
Parecem todas iguais, mas observe na mesma escala (desconsiderando as unidades de medida):
Para padronizar qualquer distribuição, basta aplicar a todos os seus valores x:
\(z = { {x - \mu} \over \sigma }\)
Subtrair a média faz com que a distribuição fique centrada em \(0\) e dividir por \(\sigma\) faz com que o desvio-padrão seja igual a \(1\). A distribuição resultante é dada em escore \(z\).
Caso normalizemos todas as curvas acima, obteremos:
A vantagem é que conhecemos as propriedades de qualquer distribuição normal, mas da normal padronizada memorizamos facilmente seus principais valores.
|
A função R que (dado o valor q, devolve a probabilidade) é: pnorm(q, mean = 0, sd = 1, lower.tail = TRUE) que calcula (por default) as probabilidades acumuladas de -infinito (lower.tail=TRUE) ao valor q solicitado. Assume, também por default, média igual a zero (mean=0) e desvio-padrão (em ingês, standard deviation) igual a 1 (sd=1), portanto uma normal padronizada. Não é necessário converter tudo para a normal padronizada quando quiser encontrar quaisquer as probabilidades. No exemplo, com média igual a 15 e desvio-padrão igual a 8, encontramos as áreas, respectivamente, para \(\pm 1 sd\):
ou, na versão da normal padronizada (que serve para simplificar):
a distribuição normal é simétrica, então achamos a probabilidade acumulada de \(-\infty\) até 1 desvio-padrão, subtraímos a metade esquerda da distribuição (encontrando a área entre 0 e 1 desvio-padrão) e, então, multiplicamos por 2. Similarmente, para \(\pm 2 sd\):
e para \(\pm 3 sd\):
Note que, entre \(\pm 2 dp\), não temos exatamente 95% da área. Para achar o valor correto, a função R que faz o reverso de pnorm (dada a probabilidade, devolve o valor q) é: qnorm(p, mean = 0, sd = 1, lower.tail = TRUE) então obtemos (deixando \(2.5\%\) em cada cauda):
e
Então, 95% da área fica, aproximadamente, no intervalo dado por \(\pm 1.96 dp\).
|
Para mulheres o valor de referência da creatinina sérica é de 0.5 a 1.1 mg/dl. Assumiremos, então, distribuição aproximadamente normal com média de 0.8 mg/dl e desvio-padrão de 0.15 mg/dl.
O seguinte código usa a função rnorm() para criar, usando um gerador de números pseudo-aleatórios (random), uma população com 10000 indivíduos:
# numero de individuos
n <- 10000
# assumindo que a distribuicao eh simetrica%
mu <- 0.8 # (0.5+1.1)/2
# assumindo que deram o intervalo de 95%
dp <- 0.15 # (1.1-0.5)/4
# cria a populacao e exibe o grafico
v <- rnorm(n,mu,dp)
d <- density(v)
plot (d, type="l",
main="Distribuicao ficticia",
xlim=c(mu-4*dp,mu+4*dp),
xlab="Creatinina",
ylab="densidade",
lty=1, lwd=2)e o seguinte código exibe a distribuição padronizada
v2 <- (v-mu)/dp
d2<- density(v2)
plot (d2, type="l",
main="Distribuicao ficticia padronizada",
xlim=c(-4, 4),
xlab="z",
ylab="densidade",
lty=1, lwd=2)Da mesma forma, assumindo-se que o valor de referência para hemoglobina em mulheres é de \(13.5 \pm 1.5 g/dl\), correspondendo estes valores à média \(\pm\) 2 desvio-padrão, e que a distribuição da hemoglobina segue aproximadamente uma distribuição normal, criamos:
# numero de individuos
n <- 10000
# media
mu <- 13.5
# desvio-padrao
dp <-1.5
# cria a populacao e exibe o grafico
v <- rnorm(n,mu,dp)
d <- density(v)
plot (d, type="l",
main="Distribuicao ficticia",
xlim=c(mu-4*dp,mu+4*dp),
xlab="Hemoglobina",
ylab="densidade",
lty=1, lwd=2) a qual é igualmente padronizada:
v2 <- (v-mu)/dp
d2<- density(v2)
plot (d2, type="l",
main="Distribuicao ficticia padronizada",
xlim=c(-4, 4),
xlab="z",
ylab="densidade",
lty=1, lwd=2)O Teorema Central do Limite (TCL) é uma das descobertas mais poderosas para as análises estatísticas. Para experimentar com ele, vamos criar uma população com uma variável fictícia que não tem distribuição normal:
# Cria populacao ficticia
source ("friendlycolor.R")
# Uma variavel qualquer
N <- 1000000 # tamanho da populacao
mu1 <- 80 # media de 2/3 da populacao
sigma1 <- 30 # desvio padrao de 2/3 da populacao
mu2 <- 180 # media de 1/3 da populacao
sigma2 <- 50 # desvio padrao de 2/3 da populacao
# set.seed(123)
# criando uma populacao ficticia
pop_valores <- round(rnorm(2*N/3, mean=mu1, sd=sigma1),0)
pop_valores <- c(pop_valores,
round(rnorm(1*N/3, mean=mu2, sd=sigma2),0))
mean_pop <- mean(pop_valores)
sd_pop <- sd(pop_valores)
# distribuicao dos valores nesta populacao ficticia
dpop_valores <- density(pop_valores)
plot (dpop_valores, main="Populacao ficticia",
xlab = "Valor", ylab = "densidade",
xlim = c(min(mu1,mu2)-4*max(sigma1,sigma2),
max(mu1,mu2)+4*max(sigma1,sigma2)),
ylim = c(0,max(dpop_valores$y)*1.3),
col = friendlycolor(10),
lwd=2, type = "l")
tp <- 44
for (i in -3:3)
{
lines(c(mean_pop-i*sd_pop, mean_pop),
c(max(dpop_valores$y)*1.2,max(dpop_valores$y)*1.2),
lwd=10, lty=1, col = paste(friendlycolor(8),tp,sep="") )
}
legend("right",
c("Distribuicao", "Media +-3 d.p."),
col=c(friendlycolor(8),paste(friendlycolor(8),"44",sep="")),
lwd=c(2,10),
lty=c(1,1),
box.lwd=0, bg="transparent") Retiraremos desta população B amostras com n indivíduos. As amostras aparecem em laranja.
B <- 3000
n <- 36
# reapresenta a populacao
plot (dpop_valores, main=paste("Amostragem (",B," amostras com n = ",n,")",sep=""),
xlab = "Valor", ylab = "densidade",
xlim = c(min(mu1,mu2)-4*max(sigma1,sigma2),
max(mu1,mu2)+4*max(sigma1,sigma2)),
ylim = c(0,max(dpop_valores$y)*1.3),
col = friendlycolor(10),
lwd=2, type = "l")
# plota media e dp populacional
tp <- 44
for (i in -3:3)
{
lines(c(mean_pop-i*sd_pop, mean_pop),
c(max(dpop_valores$y)*1.2,max(dpop_valores$y)*1.2),
lwd=10, lty=1, col = paste(friendlycolor(8),tp,sep="") )
}
amo_med <- c() # guardando as medias amostras
amo_sd <- c() # guardando os d.p. amostras
for (a in 1:B)
{
amostra <- sample(pop_valores, n, replace=FALSE)
amo_med <- c(amo_med,mean(amostra))
amo_sd <- c(amo_sd,sd(amostra))
amo_dens <- density(amostra, bw = 4)
lines(amo_dens, col=paste(friendlycolor(22),"04",sep=""), lwd=0.4)
}
mean_amo <- mean(amo_med)
sd_amo <- mean(amo_sd)
abline(v=mean_amo, lwd=2, lty=3, col=friendlycolor(19))
tp <- 44
for (i in -3:3)
{
lines(c(mean_pop-i*sd_pop, mean_pop),
c(max(dpop_valores$y)*1.2,max(dpop_valores$y)*1.2),
lwd=10, lty=1, col = paste(friendlycolor(8),tp,sep="") )
}
tp <- 44
for (i in -3:3)
{
lines(c(mean_amo-i*sd_amo, mean_amo),
c(max(dpop_valores$y)*1.3,max(dpop_valores$y)*1.3),
lwd=10, lty=1, col = paste(friendlycolor(19),tp,sep="") )
}
legend(x=mean_pop+2*sd_pop, y=max(dpop_valores$y)*1.1,
c("Distribuicao pop.",
"Media pop +-3 d.p.",
paste(B,"amostras"),
"Media(sampling)",
"Sampling +-3 d.p."
),
col=c(
friendlycolor(8),
paste(friendlycolor(8),"44",sep=""),
friendlycolor(22),
friendlycolor(19),
paste(friendlycolor(19),"44",sep="")
),
lwd=c(2,10,2,2,10),
lty=c(1,1,1,3,1),
box.lwd=0, bg="transparent",
cex=0.8) A linha pontilhada vertical em vermelho corresponde à media das medias amostrais. Na parte alta do grafico a barra horizontal em vermelho mostra a média das médias amostrais e a média dos desvios-padrão amostrais.
Repare que a média das médias amostrais coincide com a média populacional e que a média dos desvios-padrao amostrais coincide com o d.p. populacional.
v <- ""
v <- paste(v,"Populacao:\n")
v <- paste(v,"\tmedia populacional:",round(mean_pop,3),"\n")
v <- paste(v,"\td.p. populacional:",round(sd_pop,3),"\n")
v <- paste(v,"\n")
v <- paste(v,"Amostras:",B,"com n =",n,"\n")
v <- paste(v,"\tmedia das medias amostrais:",round(mean_amo,3),"\n")
v <- paste(v,"\tmedia dos d.p. amostrais:",round(sd_amo,3))
cat(v)## Populacao:
## media populacional: 113.228
## d.p. populacional: 60.41
##
## Amostras: 3000 com n = 36
## media das medias amostrais: 113.279
## media dos d.p. amostrais: 60.062
O próximo gráfico mostra a distribuição das médias amostrais.
damo_med <- density(amo_med)
plot (damo_med, main=paste("Distribuicao de Medias Amostrais\n(",B," amostras com n = ",n,")",sep=""),
xlab = "Valor", ylab = "densidade",
ylim = c(0, max(damo_med$y)*1.3),
col = friendlycolor(19),
lwd=3, type = "l")
sd_epmamostral <- sd(amo_med)
minx <- mean_amo-5*sd_epmamostral
maxx <- mean_amo+5*sd_epmamostral
byx <- (maxx-minx)/100
x_normal <- seq(from=minx, to=maxx, by=byx)
y_normal <- dnorm(x_normal, mean=mean_amo, sd=sd_epmamostral)
lines(x_normal,y_normal, lwd=2, lty=3, col = friendlycolor(19))
legend("topright",
c("Distribuicao",
"Normal"
),
col=c(
friendlycolor(19),
friendlycolor(19)
),
lwd=c(3,2),
lty=c(1,3),
box.lwd=0, bg="transparent",
cex=0.8) O TCL estabelece que a distribuição das médias amostrais, independentemente da forma da distribuição da população original tem distribuicao NORMAL com média igual à média da população e com certo desvio-padrão.
A linha pontilhada é uma distribuição normal sobreposta para mostrar que a distribuição simulada (linha sólida) é aproximadamente normal apesar da distribuição da variável original não ter este tipo de distribuição.
O desvio-padrão das médias amostrais recebe o nome de Erro Padrão das Médias Amostrais (EPM), o qual (em teoria) é estimado pelo desvio-padrão da população dividido pela raiz quadrada do tamanho (n) de cada amostra:
\(EPM = {\sigma \over \sqrt{n}}\)
|
Você pode executar o Rscript Populacao_e_EPM.R Observe o comportamento do EPM quando aumenta o tamanho da amostra. |
Como as amostras tiveram tamanho de 36, o EPM é 6 vezes menor que o desvio-padrão populacional (ou o desvio-padrão médio das 3000 amostras)
v <- ""
v <- paste(v,"\tmedia dos d.p. amostrais:",round(sd_amo,3),"\n")
v <- paste(v,"\td.p. das médias amostrais (EPM):",round(sd_epmamostral,3),"\n")
cat(v)## media dos d.p. amostrais: 60.062
## d.p. das médias amostrais (EPM): 9.985
Ninguem faz B amostras de uma população. Na prática somente uma amostra é obtida.
# reapresenta a populacao
plot (dpop_valores, main=paste("Populacao e amostra unica com n = ",n,sep=""),
xlab = "Valor", ylab = "densidade",
xlim = c(min(mu1,mu2)-4*max(sigma1,sigma2),
max(mu1,mu2)+4*max(sigma1,sigma2)),
ylim = c(0,max(dpop_valores$y)*1.3),
col = friendlycolor(10),
lwd=2, type = "l")
# amostra unica
amostra_unica <- sample(pop_valores, n, replace=FALSE)
mean_amouni <- mean(amostra_unica)
sd_amouni <- sd(amostra_unica)
amouni_dens <- density(amostra_unica, bw = 4)
lines(amouni_dens, col=paste(friendlycolor(22),sep=""), lwd=3)
abline(v=mean_amouni, lwd=2, lty=3, col=friendlycolor(19))
tp <- 44
for (i in -3:3)
{
lines(c(mean_pop-i*sd_pop, mean_pop),
c(max(dpop_valores$y)*1.2,max(dpop_valores$y)*1.2),
lwd=10, lty=1, col = paste(friendlycolor(8),tp,sep="") )
}
tp <- 44
for (i in -3:3)
{
lines(c(mean_amouni-i*sd_amouni, mean_amouni),
c(max(dpop_valores$y)*1.3,max(dpop_valores$y)*1.3),
lwd=10, lty=1, col = paste(friendlycolor(19),tp,sep="") )
}
legend(x=mean_pop+2*sd_pop, y=max(dpop_valores$y)*1.1,
c("Distribuicao pop.",
"Media pop +-3 d.p.",
"Amostra unica",
"Media amostral",
"Media amostral +-3 d.p."
),
col=c(
friendlycolor(8),
paste(friendlycolor(8),"44",sep=""),
friendlycolor(22),
friendlycolor(19),
paste(friendlycolor(19),"44",sep="")
),
lwd=c(2,10,2,2,10),
lty=c(1,1,1,3,1),
box.lwd=0, bg="transparent",
cex=0.8) A amostra, se não houver problemas, traz propriedades da população:
v <- ""
v <- paste(v,"Populacao:\n")
v <- paste(v,"\tmedia populacional:",round(mean_pop,3),"\n")
v <- paste(v,"\td.p. populacional:",round(sd_pop,3),"\n")
v <- paste(v,"\n")
v <- paste(v,"Amostra com n =",n,"\n")
v <- paste(v,"\tmedia amostral:",round(mean_amouni,3),"\n")
v <- paste(v,"\td.p. amostral:",round(sd_amouni,3),"\n")
cat(v)## Populacao:
## media populacional: 113.228
## d.p. populacional: 60.41
##
## Amostra com n = 36
## media amostral: 103.139
## d.p. amostral: 58.833
Na prática, porém, não temos a população como referência e, assim, não fazemos ideia se temos uma amostra representativa:
# reapresenta sem a populacao
plot (NA, main=paste("Amostra unica com n = ",n,sep=""),
xlab = "Valor", ylab = "densidade",
xlim = c(min(mu1,mu2)-4*max(sigma1,sigma2),
max(mu1,mu2)+4*max(sigma1,sigma2)),
ylim = c(0,max(dpop_valores$y)*1.3)
)
# reapresenta a ultima amostra
lines(amouni_dens, col=paste(friendlycolor(22),sep=""), lwd=3)
abline(v=mean_amouni, lwd=2, lty=3, col=friendlycolor(19))
tp <- 44
for (i in -3:3)
{
lines(c(mean_amouni-i*sd_amouni, mean_amouni),
c(max(dpop_valores$y)*1.3,max(dpop_valores$y)*1.3),
lwd=10, lty=1, col = paste(friendlycolor(19),tp,sep="") )
}
legend(x=mean_pop+2*sd_pop, y=max(dpop_valores$y)*1.1,
c("Amostra unica",
"Media amostral",
"Media amostral +-3 d.p."
),
col=c(
friendlycolor(22),
friendlycolor(19),
paste(friendlycolor(19),"44",sep="")
),
lwd=c(2,2,10),
lty=c(1,3,1),
box.lwd=0, bg="transparent",
cex=0.8) O que podemos fazer é um processo similar ao anterior, fazendo B reamostragens de n elementos. Isto parece inútil, pois sempre teremos a própria amostra:
amostra <- c(1, 2, 3, 4, 5)
n <- length(amostra)
cat("amostra: ", amostra,
": média =", round(mean(amostra),2),
"e dp =", round(sd(amostra),2),
"\n")## amostra: 1 2 3 4 5 : média = 3 e dp = 1.58
# 9 reamostragens sem reposicao
for (r in 1:9)
{
# *** replace=FALSE por default ***
reamostra <- sample(amostra,n)
cat("reamostra",r,": ", reamostra,
": média =", round(mean(reamostra),2),
"e dp =", round(sd(reamostra),2),
"\n")
}## reamostra 1 : 2 5 4 3 1 : média = 3 e dp = 1.58
## reamostra 2 : 5 2 4 1 3 : média = 3 e dp = 1.58
## reamostra 3 : 5 4 3 1 2 : média = 3 e dp = 1.58
## reamostra 4 : 2 5 1 3 4 : média = 3 e dp = 1.58
## reamostra 5 : 3 4 1 5 2 : média = 3 e dp = 1.58
## reamostra 6 : 2 5 1 4 3 : média = 3 e dp = 1.58
## reamostra 7 : 2 5 1 4 3 : média = 3 e dp = 1.58
## reamostra 8 : 5 1 3 2 4 : média = 3 e dp = 1.58
## reamostra 9 : 5 3 4 1 2 : média = 3 e dp = 1.58
Então, para obtermos variações de nossa única amostra, fazemos a reamostragem COM REPOSICAO:
# 9 reamostragens com reposicao
for (r in 1:9)
{
# *** replace=TRUE permite a reposicao ***
reamostra <- sample(amostra,n,replace=TRUE)
cat("reamostra",r,": ", reamostra,
": média =", round(mean(reamostra),2),
"e dp =", round(sd(reamostra),2),
"\n")
}## reamostra 1 : 5 5 4 5 2 : média = 4.2 e dp = 1.3
## reamostra 2 : 2 1 5 5 3 : média = 3.2 e dp = 1.79
## reamostra 3 : 2 4 3 4 3 : média = 3.2 e dp = 0.84
## reamostra 4 : 3 2 5 1 4 : média = 3 e dp = 1.58
## reamostra 5 : 1 5 2 2 5 : média = 3 e dp = 1.87
## reamostra 6 : 5 4 4 1 1 : média = 3 e dp = 1.87
## reamostra 7 : 2 5 1 2 1 : média = 2.2 e dp = 1.64
## reamostra 8 : 2 2 3 3 2 : média = 2.4 e dp = 0.55
## reamostra 9 : 2 2 2 5 2 : média = 2.6 e dp = 1.34
Voltando à população, podemos fazer as B amostras com n elementos com reposição:
B <- 3000
n <- 36
# reapresenta a amostra
plot (NA, main=paste("Bootstrapping (",B," reamostras com n = ",n,")",sep=""),
xlab = "Valor", ylab = "densidade",
xlim = c(min(mu1,mu2)-4*max(sigma1,sigma2),
max(mu1,mu2)+4*max(sigma1,sigma2)),
ylim = c(0,max(dpop_valores$y)*1.3),
col = friendlycolor(22),
lwd=2, type = "l")
boot_med <- c() # guardando as medias amostras do bootstrapping
boot_sd <- c() # guardando os d.p. amostras do bootstrapping
for (a in 1:B)
{
amostra <- sample(amostra_unica, n, replace=TRUE)
boot_med <- c(boot_med,mean(amostra))
boot_sd <- c(boot_sd,sd(amostra))
amo_dens <- density(amostra, bw = 4)
lines(amo_dens, col=paste(friendlycolor(15),"04",sep=""), lwd=0.4)
}
# replota a amostra unica
lines (amouni_dens, lwd=3, col = friendlycolor(22))
mean_boot <- mean(boot_med)
sd_boot <- mean(boot_sd)
abline(v=mean_boot, lwd=2, lty=3, col=friendlycolor(13))
# replota media e dp da amostra unica
tp <- 44
for (i in -3:3)
{
lines(c(mean_amouni-i*sd_amouni, mean_amouni),
c(max(dpop_valores$y)*1.3,max(dpop_valores$y)*1.3),
lwd=10, lty=1, col = paste(friendlycolor(19),tp,sep="") )
}
# plota media e dp do bootstrapping
tp <- 44
for (i in -3:3)
{
lines(c(mean_boot-i*sd_boot, mean_boot),
c(max(dpop_valores$y)*1.2,max(dpop_valores$y)*1.2),
lwd=10, lty=1, col = paste(friendlycolor(13),tp,sep="") )
}
legend(x=mean_pop+2*sd_pop, y=max(dpop_valores$y)*1.1,
c("Dist. da amostra",
"Media am +-3 d.p.",
paste(B,"amostras"),
"Media geral (bootstrp.)",
"Medias bootstrp",
"Bootstrp. +-3 d.p."
),
col=c(
friendlycolor(22),
paste(friendlycolor(19),"44",sep=""),
friendlycolor(15),
friendlycolor(13),
paste(friendlycolor(24),"20",sep=""),
paste(friendlycolor(13),"44",sep="")
),
lwd=c(2,10,2,2,10,10),
lty=c(1,1,1,3,1,1),
box.lwd=0, bg="transparent",
cex=0.8) E verificar que, novamente, a distribuição das médias amostrais é normal:
# apresenta o EPM obtido por bootstrapping
dboot_med <- density(boot_med)
plot (dboot_med, main=paste("Distribuicao de Medias Amostrais\n(",B," amostras com n = ",n,")",sep=""),
xlab = "Valor", ylab = "densidade",
ylim = c(0, max(damo_med$y)*1.3),
col = friendlycolor(15),
lwd=3, type = "l"
)
abline(v=mean_boot, lwd=2, lty=3, col=friendlycolor(15))
sd_epmboot <- sd(boot_med)
minx <- mean_boot-5*sd_epmboot
maxx <- mean_boot+5*sd_epmboot
byx <- (maxx-minx)/100
x_normal <- seq(from=minx, to=maxx, by=byx)
y_normal <- dnorm(x_normal, mean=mean_boot, sd=sd_epmboot)
lines(x_normal,y_normal, lwd=2, lty=3, col = friendlycolor(15))
icx <- quantile (boot_med, probs=c(0.025,0.975))
icy <- rep(max(damo_med$y)*1.3,3)
icy <- icy * c(0.95,0.9,0.85)
lines(c(icx[1],icx[1],icx[1],icx[2],icx[2],icx[2]),
c(icy[1],icy[3],icy[2],icy[2],icy[1],icy[3]),
lwd=3, col = friendlycolor(15))
points(mean_boot, icy[2], pch=21, bg = friendlycolor(45), col = friendlycolor(15), cex=2)
points(mean_pop, icy[2], pch=21, bg = friendlycolor(19), col = friendlycolor(19), cex=2)
legend("topright",
c("Dist. bootstrapping",
"Normal bootstrapping",
"Media populacao",
"Media bootstrapping"
),
col=c(
friendlycolor(15),
friendlycolor(15),
friendlycolor(19),
friendlycolor(15)
),
pt.bg=c(NA,NA,friendlycolor(19),friendlycolor(45)),
lwd=c(3,2,NA,NA),
lty=c(1,3,NA,NA),
pch=c(NA,NA,21,21),
box.lwd=0, bg="transparent",
cex=0.8) O valor de EPM, a partir do bootstrapping, também é próximo a 1/6 do desvio-padrão populacional.
v <- ""
v <- paste(v,"Populacao (que nunca veremos):\n")
v <- paste(v,"\tmedia populacional:",round(mean_pop,3),"\n")
v <- paste(v,"\td.p. populacional:",round(sd_pop,3),"\n")
v <- paste(v,"\n")
v <- paste(v,"Amostra com n =",n,"\n")
v <- paste(v,"\tmedia amostral:",round(mean_amouni,3),"\n")
v <- paste(v,"\td.p. amostral:",round(sd_amouni,3),"\n")
v <- paste(v,"Reamostras:",B,"com n =",n,"\n")
v <- paste(v,"\tmedia das medias amostrais:",round(mean_boot,3),"\n")
v <- paste(v,"\tmedia dos d.p. amostrais:",round(sd_boot,3),"\n")
v <- paste(v,"EPM por bootstrapping a partir da amostra:","\n")
icx <- round(icx,3)
v <- paste(v,"\tmedia das medias amostrais:",round(mean_boot,3),"[",icx[1],",",icx[2],"]","\n")
v <- paste(v,"\tEPM (desvio-padrao das medias amostrais):",round(sd_epmboot,3),"\n")
cat(v)## Populacao (que nunca veremos):
## media populacional: 113.228
## d.p. populacional: 60.41
##
## Amostra com n = 36
## media amostral: 103.139
## d.p. amostral: 58.833
## Reamostras: 3000 com n = 36
## media das medias amostrais: 103.121
## media dos d.p. amostrais: 57.198
## EPM por bootstrapping a partir da amostra:
## media das medias amostrais: 103.121 [ 85.388 , 123.056 ]
## EPM (desvio-padrao das medias amostrais): 9.569
Observe, ainda, que a média populacional está dentro do intervalo de confiança de 95% estimado pelo bootstrapping.